home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Ham Radio 2000
/
Ham Radio 2000.iso
/
ham2000
/
misc
/
radiosky
/
azel.bas
next >
Wrap
BASIC Source File
|
1996-05-15
|
5KB
|
126 lines
1 'The following simple program, written in IBM BASIC, will perform the
2 'calculations to convert HA, decl to and from az-el, including the
3 'computation of LST. Although this program prints out the results
4 'to limited precision, the calculations are precise to a small
5 'fraction of a degree, and are valid at least throughout this and
6 'the next century. SKYCOORD.BAS is a complete working program, but
7 'the reader may wish to use it as a template for his or her own program
8 'in any favorite programming language, perhaps adding a more sophisticated
9 'user interface.
10 REM Convert RA, Dec to Az, El and vice versa.
20 REM Put together by Darrel Emerson, aa7fv, July 1995.
25 REM From various sources, especially "Astronomy with your PC"
26 REM by Peter Duffet-Smith, Cambridge University Press.
30 DEFDBL A-Z
40 pi = 4# * ATN(1#)
50 INPUT "Day, Month, Year (e.g. 17,6,1995)"; DY, MN, yr
70 INPUT "Observer's longitude (decimal degrees, W is negative)"; lo
80 gl = lo * pi / 180#
90 INPUT "Observer's latitude (decimal degrees, S is negative)"; la
100 tl = la * pi / 180#
120 PRINT
130 PRINT "RA-Dec to Az-El (1)": PRINT "Az-El to RA-Dec (2)"
140 INPUT " or finish. (3) "; ir
150 IF ir < 1 OR ir > 2 THEN PRINT "End conversion program": STOP
160 PRINT : INPUT "UT hours (decimal hours: e.g. 23:45 = 23.75)"; ut
170 GOSUB 2100: REM calculate LST hours
180 PRINT "LST"; USING "####.####"; sl; : PRINT " hours"
190 IF ir = 1 THEN GOSUB 300
200 IF ir = 2 THEN GOSUB 500
210 GOTO 120
300 REM Convert RA-Dec to Az-El
310 INPUT "Right Ascension, decimal hours (e.g. 20h54m=20.9)"; rh
320 INPUT "Declination, decimal degrees"; dd
330 de = dd * pi / 180#
340 ha = (sl - rh) * 15# * pi / 180#
350 GOSUB 3000
360 PRINT "Azimuth: "; USING "####.#"; az * 180# / pi; : PRINT " degrees"
370 PRINT "Elevation:"; USING "####.#"; el * 180# / pi; : PRINT " degrees"
380 RETURN
500 REM Convert Az-El to RA-Dec
510 INPUT "Azimuth, decimal degrees"; ad
520 INPUT "Elevation, decimal degrees"; ed
530 az = ad * pi / 180#
540 el = ed * pi / 180#
550 GOSUB 4000
560 hh = ha * 180# / pi / 15#: REM convert hour angle into hours
570 REM print"Hour Angle:";using "###.####";hh;:print " hours"
580 ra = sl - hh
590 IF ra < 0# THEN ra = ra + 24#
600 IF ra > 24# THEN ra = ra - 24#
610 PRINT "Right Ascension: "; USING "###.##"; ra; : PRINT " hours"
620 PRINT "Declination: "; USING "###.#"; de * 180 / pi; : PRINT " degrees"
630 RETURN
1000 REM subroutine atan2, using sine s and cosine c.
1005 REM All angles radians.
1010 IF c = 0 AND s = 0 THEN atan2 = 0: RETURN
1020 IF c = 0 AND s > 0 THEN atan2 = pi / 2#: RETURN
1030 IF c = 0 AND s < 0 THEN atan2 = -pi / 2#: RETURN
1040 atan2 = ATN(s / c): IF c < 0 THEN atan2 = atan2 + pi
1050 IF atan2 > pi THEN atan2 = atan2 - 2# * pi
1060 REM atan2 is returned between +/-pi
1070 RETURN
1100 REM subroutine to calculate Julian days since the start of the
1105 REM century. (i.e. since 1900 Jan 0.5)
1110 REM TAKES YR (E.G. 1995), MN (E.G. 6), DY (E.G. 17) & gives DJ
1120 M1 = MN: Y1 = yr: B0 = 0
1130 IF MN < 3 THEN M1 = MN + 12: Y1 = Y1 - 1
1140 A0 = INT(Y1 / 100): B0 = 2 - A0 + INT(A0 / 4)
1150 C0 = INT(365.25# * Y1) - 694025!
1160 D0 = INT(30.6001# * (M1 + 1)): dj = B0 + C0 + D0 + DY - .5#
1170 RETURN
2100 REM Find GST (=sg) and LST (=sl) from UT,
2105 REM using UT, YR, MN, DY, GL (=longitude, W neg)
2110 GOSUB 1100: REM get DJ, days since 1900
2120 d = INT(dj - .5#) + .5#: t = (d / 36525#) - 1
2130 r1 = 6.697374558# + (2400# * (t - ((yr - 2000#) / 100#)))
2140 r0 = t * (.0513366# + t * (.00002586222# - t * .000000001722#))
2150 t0 = (r0 + r1) - 24# * INT((r0 + r1) / 24#)
2160 sg = (ut * 1.002737908#) + t0
2170 sg = sg - 24# * INT(sg / 24#)
2180 sl = sg + (gl * 180# / pi / 15#):
2190 sl = sl - 24# * INT(sl / 24#)
2200 REM sg is Grenwich Sidereal Time, sl is Local Sidereal Time
2210 RETURN
3000 REM Find Azimuth AZ and Elevation EL.
3010 REM from Hour Angle HA, Declination DE and Terrestrial Latitude TL.
3020 REM All angles radians.
3030 c1 = -COS(de) * COS(ha) * SIN(tl) + SIN(de) * COS(tl)
3040 c2 = -COS(de) * SIN(ha)
3050 c3 = COS(de) * COS(ha) * COS(tl) + SIN(de) * SIN(tl)
3060 s = c2: c = c1: GOSUB 1000: az = atan2
3070 s = c3
3080 IF ABS(SIN(az)) > .7 THEN c = c2 / SIN(az): GOTO 3110
3090 c = c1 / COS(az)
3100 REM The "0.7" is a criterion for minimizing numerical rounding error
3110 GOSUB 1000: el = atan2
3120 IF az < 0 THEN az = az + 2# * pi
3130 REM "az" now between 0 and 2*pi, with Az=0 north
3140 RETURN
4000 REM Find Hour Angle HA, and Declination DE,
4010 REM from Azimuth AZ and Elevation EL, with Latitude TL.
4020 REM All angles radians.
4030 c1 = -COS(el) * COS(az) * SIN(tl) + SIN(el) * COS(tl)
4040 c2 = -COS(el) * SIN(az)
4050 c3 = COS(el) * COS(az) * COS(tl) + SIN(el) * SIN(tl)
4060 s = c2: c = c1: GOSUB 1000: ha = atan2
4070 s = c3
4080 IF ABS(SIN(ha)) > .7 THEN c = c2 / SIN(ha): GOTO 4110
4090 c = c1 / COS(ha)
4100 REM The "0.7" is a criterion for minimzing numerical rounding error
4110 GOSUB 1000: de = atan2
4120 IF ha < 0 THEN ha = ha + 2# * pi
4130 REM ha now between 0 and 2*pi
4140 RETURN